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Q^ . Two-dimensional decaying turbulent flow is known to approach apparently stable 

states after a long time evolution. A few theories and models have been so far pro- 
posed to account for this relaxation. In this paper, we compare results of numerical 
experiments with the predictions of these theories to assess their applicability. We 
O I study the long time decay of initially multilevel vorticity fields on the periodic box, 

^ ' and characterize the outcoming final states. Our final states do not match the pre- 

O ■ dictions of the theories; a broader variety of dipole profiles, as well as nonstationary 

final states are found. The problem of the robustness of the relaxational state with 
k> ^ respect to variations of the Reynolds number and different numerical resolution is 

j_j ■ addressed. The observed configurations also do not necessarily possess the maximal 
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energy, in contrast to what is anticipated by some of the theories. We are led to 
conclude that the mixing of the vorticity is generally not ergodic, and that some 
metastable configurations may inhibit the attainment of an equilibrium state. 
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1 Introduction and motivations 



The decay of the vorticity field of a two-dimensional incompressible fiuid, 
which obeys to the unforced Navier-Stokes equation, shows several interest- 
ing features. Stable large scale structures may form, organize themselves on 
the scale of the accessible domain, and survive for a long time before being 
ultimately damped by the dissipation. Occurrences of such structures are ob- 
served in numerical and laboratory experiments, as well as in planetary scale 
fiows. We will call from here on such configurations "final" ; by that we mean 
subjected only to viscous decay and no more to filamentation and mixing. 
Provided that the Reynolds number of the flow is large enough, the timescales 
for these two processes can be very different. We are concerned with the ap- 
pearance of final states emerging from arbitrary initial conditions, with their 
robustness and their predictability. Most of the numerical simulations pre- 
sented in the literature start from random initial conditions, whose Fourier 
spectra decay variously in k, with uncorrelated phases. In such cases, inter- 
mediate states with few isolated vortex cores are observed. In the traditional 
scenario (Benzi et al. 1988), the progressive merging of likely signed vortex 
cores is then observed. Several papers in the literature deal with the merging 
process, either analyzing the evolution in time of the vortex census, or mod- 
elling the process with simplified dynamics. Less attention has been so far 
devoted to the actual characteristics of the latest state of the field. 

A few explanations for the occurrence of final states have been proposed so far. 

We will refer specifically to those theories, which are formulated in the physical 

space rather than in the Fourier one: the minimum enstrophy assumption and 

the maximal entropy state statistical theory (Miller ct al. 1992, Robert and Sommeria 1991). 

Some numerical evidence has also been invoked to reconnect computed results 

with the Joyce- Montgomery equation, which is originally derived for the mean 

field hmit of a system of point vortices (Montgomery et al. 1993). 

We consider several cases of decay of vorticity fields. All cases start from initial 
conditions which should be generic and appropriate for the application of the 
final state theories. Our results contrast with the conclusion that an universal 
state emerges out of the relaxation, without memory of the dynamical path 
which led to it. In fact, we even see that the final 'state' may be unsteady: 
some fields attain a late time configuration which is not stationary. Such late 
fields may move quasiperiodically and steadily, being just slowly damped by 
the (arbitrary small) dissipation. 

This paper is organized as follows: in section 2 we briefly review the treatment 
of the problem of two-dimensional turbulent decay, highlighting the different 
variational problems which are solved to find the final state, and the functional 
relations ci;(^/') derived. In section 3 we describe the numerical experiments per- 
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formed. Finally, in section 4 we comment on the open aspects which require 
further investigation. A number of remarks which are related to the maxi- 
mization of the energy and its relevance for the final state are left for the 
Appendix. 



2 Review of the available theories 

Let us first recall the notation. As is well known, the equation of motion for 
the vorticity ci;(x) is written as 



u;t{x) = -v(x) • Vu;(x) + i/V^a;(x) = J{uj{x),iJj{x)) + uW^u;{x) . (1) 



In (1), i/j denotes the streamfunction, obtained from a;(x) by means of the 
Green function G of the Laplacian operator 



The integral is carried on the fiuid domain with the proper boundary condi- 
tions, so that V^V = Introducing the notation V"*" = {dy, —dx), we may 
write V-^ip = v and J(a;, ip) = Va; ■ V-^ip. The energy of the field is defined as 



with / a positive integer. The quantity Q2{= Q) is traditionally called the 
enstrophy. Another global quantity which is defined is the palinstrophy 




(2) 




(3) 



and the moments of the vorticity as 




(4) 




(5) 



These quantities evolve in time according to: 



Et = -2vQ , Qt = -2vP , 
Qi,t = - l)v I u;(x)'-2 [Va;(x)]' d'x 



(6) 
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E and all Qi are constants of motion if z/ = 0. In the limit of vanishing 
viscosity, is a constant of motion, but Q may not be, because P can get 
larger inversely proportional to i^. 

We refer to theories which predict the final state as the most probable outcome 
of the decay. In a way or the other, all these models assume a distinction 
between the fully detailed dynamics expressed by equation (1), and that of 
a reduced set of macroscopically observable quantities. The evolution of the 
field is seen as a process in which the initial information is lost in some way. 
As a deliberate simplification, the final state is sought as the one which is 
fully described only by a few macroscopical constraints, which are often called 
"rugged invariants" . Such theories treat the viscous, finite resolution problem, 
as one in which some quantities are "better" conserved than the rest, in place 
of the infinite set of the inviscid case. There is indeed some ambiguity, and 
properly speaking the case = is different from the limit — > 0, since in 
the former infinitely steep gradients might form. Where possible, methods of 
statistical mechanics are applied, and some extremum principle is invoked. A 
comprehensive review of the various positions can be found elsewhere (e.g. 
Miller et al. 1992); here we recall them briefiy, and discuss their conclusions. 

2. 1 Equilibrium Fourier spectra 

A first group of theories is formulated in the Fourier space. The older Krai- 
chnan-Batchelor-Leith statistical theory (Kraichnan 1967) predicts a Fourier 
spectrum E{k) ~ k~^, relying on the assumption of the locality of the interac- 
tions among the components. An improvement by Kraichnan (1971), based on 
the test-field-model closure approximation, corrected this spectrum to E{k) ~ 

(In /c//co)^- While some numerical simulations (Kida et al. 1988,Brachet et al. 1988) 
support these spectra, quite different spectra have been observed by others 
(see for example McWilliams 1984) The reason of the discrepancy is not clear, 
though the formation of stable structures may take a key role. 

A later theory due to Kraichnan (Kraichnan and Montgomery 1980,Carnevale 1982) 
proposes a statistical mechanics for the energies of the Fourier components. 
An ultaviolet cutoff in k has to be enforced. Only E and Q are assumed to be 
constants of motion, and are fixed as constraints. No correlation is assumed 
between the phases of a; (A;), and no other moment of the vorticity is conserved. 
This theory predicts a statistical equilibrium spectrum 

with arbitrary constants a and (5. The agreement of these spectra with those 
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coming from numerical simulations, and especially with ours, is controversial. 



2.2 Point vortex systems 

A second line of reasoning considers the statistical properties of an ensemble of 

point vortices. The rationale for connecting vorticity fields with such systems 

is that a system of point vortices approximates weakly, in the continuum limit, 

the Euler equation (Caglioti et al. 1992,Eyink and Spohn 1993, Campbell and O'Neil 1991); 

the full Navier-Stokes equation can be emulated by point vortices which diffuse 

with an additional Brownian motion (Chorin 1994). An entropy of the system 

is introduced and maximized. In the mean field limit, a differential equation is 

derived for the equilibrium configuration (Montgomery and Joyce 1974,Kida 1975,Pointin and Lundgre 



This provides us with a first example of an equation which relates functionally 
ujo and ipo. As is known, the functional dependence implies the stationarity of 
the motion, in the case of null dissipation. In the special case of an equal num- 
ber of opposite charged positive and negative vortices, the Joyce-Montgomery 
equation reduces to the sinh-Poisson equation 



This equation has been furthermore studied referring to the inverse scattering 
theory for the sin-Gordon equation (Ting et al. 1987). Solutions possessing 
simple scattering spectra can be constructed, but no dynamical analysis, be- 
yond a comparison of shapes, was done. 

Montgomery et al. (1993) give an interpretation of the problem which is some- 
how related. They propose a decomposition of the vorticity field in four non- 
physical positive subfields. An entropy of the form 



is then maximized individually for each subfield, with the proper constraints 
on the total energy and vorticity. The remaining arbitrary constants are fitted 
to the results of a single high resolution, high Reynolds number Navier-Stokes 
numerical simulation. For that case, they achieve a good fit of the ip{uj) scatter- 
plot at late times. For our purposes, it suffices to note that their final relation 
implies 






(9) 




(10) 



u;o(x) = ci e' 



C2& 



+C3, 




5 



which is an elaboration of (8), and can be assimilated to equation (21) in the 
case of 3 levels. 



2.3 Minimum enstrophy principle 

The identification of the final state as the one with lowest enstrophy dates back 

to Bretherton and Haidvogel (1976). They argued that while the energy can 

almost be conserved by a good numerical scheme, the vorticity filamentates 

progressively and smoothes out. Arguments related to the universality of the 

energy-enstrophy cascades, as in the Kraichnan-Batchelor-Leith theory, would 

predict for instance a behavior of Q(t) ~ t"'^ (Carnevale et al. 1992,Bartello and Warn 1996). 

The idea of a faster decay of the enstrophy with respect to the energy is often 

referred to as the "selective decay hypothesis" . The final state is consequently 

found variationally, by minimizing Q with constrained E. According to this 

hypothesis, axisymmetric vortex shapes on the infinite plane can be calculated 

(Leith 1984). Time asymptotic estimates for closed square box solutions are 

discussed by van Groesen (1988). 

In the case of doubly periodical boundary conditions, it is straightforward to 
find a solution. Imposing 



which on the periodic square admits solutions of the form ujo{x.) = X^^^ke''^'^, 
with |kp = A equal to a squared integer. This implies that Q = XE. For a 
given energy the minimal enstrophy is then achieved for A = 1, and the general 
sohition becomes a;o(x) = a;icos(x + a) + uj2Cos{y + b), with ui, lj2, a and 
b being arbitrary constants. Linear combinations of such sinusoidal solutions 
with different k (i.e. complete Fourier series) do not satisfy the requirement. 
In the numerical experiments we find final states of completely different forms. 
We remark however that this principle was introduced for flows with additional 
"topographic" terms in (1), and we do not exclude that it may provide realistic 
results in cases where these are dominant. 




(12) 



we obtain 



^^o(V'o) = Xipo , 



(13) 
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2-4 Vortex censuses and punctuated dynamics 

In a number of papers, appeared around 1990, the late lowering of the enstro- 

phy is solely explained as a result of progressive vortex mergings (McWiUiams 1990,Matthaeus et al. 19^ 
These papers consider situations with an intermediate time dynamics domi- 
nated by many well separated vortex cores, which behave approximately like 
point vortices. The subsequent evolution is schematized by a progressive col- 
lapse and merging of these objects. Statistics of the number N of cores in time, 
models for the probability of merging are examined, and eventually lead to dif- 
ferent scalings of Q{t). The final state is assumed to be a dipole, and its prop- 
erties are sought by scaling the relevant quantities down to = 2. A popular 
model is the punctuated Hamiltonian system (Carnevale et al. 1991), which 
is the traditional point vortex model fitted with a nonconservative merging as 
vortices get close enough. Such models can be further elaborated accounting 
for extended cores (Riccardi et al. 1995). 

2.5 Maximum entropy theory 

The reasoning is based on the combinatorics of infinitesimal vorticity patches, 
at a scale smaller than that which determines the energy of the configuration. 
We follow the notation of Miller et al. (1992), rather than the equivalent one 
of Robert and Sommeria (1991). The theory is indeed intended only for the 
Euler equation; Weichman (1993) proposes an additional argument in order 
to include the viscosity in connection with the underresolution, which seems 
incorrect. The theory mimics the statistical mechanics of a many particle 
system. A given distribution of (infinitesimally grained) vorticity is assimilated 
to a microstate; at the macroscopical level, only a coarse averaged vorticity 
a;(x) can be observed. The fine scale structure is remembered by introducing 
a local probability distribution n(x, (T)d(T of vorticity, which says how large is 
the probability of having a microscale vorticity cr < a;(x) < a + da at point 
X. The macroscopic averaged vorticity is then 



The macrostate is the field of probability over the whole domain. Any mi- 
croscale distribution of vorticity, which looks on the macroscale like that prob- 
ability, is said to be a compatible microstate. The fluid is expected to relax 

to the macrostate which can be achieved in the largest number of ways, and 
thus is maximally probable. The theory relies on the strong assumption, that 
the microscale mixing of vorticity is ergodic. This means that the available 
vorticity is completely free to mix in any possible (area preserving) way, so 




(14) 
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that only the probabihty determines which outcome is hkely to be observed. 
A statistical mechanics canonical approach is undertaken. A free energy 



F{{n}) = -S{{n}) - PE{{n}) + ^ i^iQi. (15) 

1=0 

is maximized. Here S{{n}) is the entropy function 

S{{n}) — — j n{yi, a) In n{x, a) (fxda , (16) 

the energy, expressed in terms of n(x, a), is 

E{{n}) = ^J aa'n{^, a)n{^', (t')G(x', ^)(fy:(P:>c'dada' , (17) 



and the constants of motion are multiplied by appropriate Lagrange factors 
and added to F({n}). 

The constraints to be imposed are E = constant, / nda = 1 and / nd^x = 
g{(j). The function g{(T) is the global vorticity distribution, which should be 
invariant for microscale inviscid flows. To implement the last constraint, the 
conservation of all moments of the vorticity is instead required. It is assumed 
to be sufficient that the infinite set of moments of the vorticity are conserved 
without requiring the topological correctness of the flux. In other words, it is 
assumed that the area preserving vorticity mappings are dense (at least in the 
coarse average) in the topologically feasible fluxes. Using g{a) instead of Qi, 

oo „ oo „ „ 

^/i.; / Lu'-d'^x = "^1^1 a^g{a)da = I iJ,{a)n{x, a)dad'^x . (18) 
1=0 1=0 



Functional derivation with respect to n(x, a) and algebraic manipulation lead 
to the system 



g-/3(o-V'o(x)-/:i(o-)) 

a;o(x) = j (7no(x, cr)d(7 = -V^Vo(x) , (19) 

which has to be solved in order to find the maximally probable macrostate 
no(x, o"). In this procedure, the dependencies ol j3 on E and of /i((T) on g{a) 
are left as implicit; their actual form is supposed to be found only after solv- 
ing consistently the system. It is also assumed, but not proven, that Lagrange 
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multipliers can be determined for any physically accessible values of the con- 
served quantities. The system (19) does not in fact say very much. For /5 < 0, 
it states that the probability of having at x a vorticity of the same sign of 
■V^o(x) grows with cr, but is shaped by the weight factor exp[/5/i.((T)]. For /? > 0, 
the same applies to a vorticity opposite in sign to ■00 (x). The fact that /3 < 
in physical situations is inferred in comparison with the case of a 2D Coulomb 
gas (Miller et al. 1992). 

Robert and Sommeria (1992) also derive an evolution equation for the ap- 
proach to the maximum entropy state in this framework. Going further on, 
Jiittner et al. (1996) propose a way to include random forcing in the maximum 
entropy theory. 

Particular forms of '4'o{u!o) can be found only with additional hypotheses. If, for 
instance, the vorticity takes only N different values {cui, . . . , cun), then n(x, a) 
will be everywhere a sum of delta functions. The system (19) becomes 



»„(x, a) = = - ■^.) ■ (20) 

ujoix) — — — -. — -— - — , — ni(-x.)uji . (21) 
^1=1 1=1 

The latter equation expresses a single valued relation a;o(^o)- Inverted as 
ipo{u!o), it may eventually be multiple-branched. Specifically, for two opposite 
levels coi — —102-, (21) becomes 

c.o(x) = tanh f-/3c.,^o(x) + f^t^A^Azl^ . (22) 



The slope of the curve in the (a;, ■0) plane is determined by /3, which depends 
on the energy; the position of the origin is fixed by //(a;i) — /i(— cui), which 
in turn can be expressed as a function of Q (all the higher moments of the 
vorticity are related to Q by the limitation to two levels). 

For a two-level vorticity distribution it is even possible to reconstruct n(x, cr) 
knowing a7(x). This happens since ni(x) and 712 (x) are found from 

ni(x) n2(x) = 1 , a;ini(x) a;2n2(x) = a;(x) . (23) 



The entropy S can thus be evaluated directly in terms of a;(x): 
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, _ J [uJ2 - (^(x)] ln[uj2 - ^(x)] + [a;(x) - cui] ln[a;(x) - cui] 



UJ2 — OJl 



+ 



(24) 



A variety of ipo{'-^o) other than (22) can be derived as well. Assuming a priori 
that PijL{a) — —\a\/q, that is, if the weight factor is Poissonian (Pasmanter 1994), 



Assuming instead that I3h{(j) = — (a/g)^, that is, if the weight factor is Gaus- 
sian (Miller et al. 1992, section VI D), 

a;o(x) = -/3gVo- (26) 



A primary importance has been attributed to the 'dressed vorticity corohary' 
(DVC) (MiUcr ct al. 1992). This corollary says that if one "guesses" from 
the macroscopical equilibrium state, that is if one writes no(x, cr) = d{a — 
a;o(x)), then a frozen dynamics is obtained, which is the (3 — > — oo limit of the 
true one. It is to remark that the "dressed vorticity distribution" 

gd{a)^ ls{a-uJo{x.))d^^ (27) 



is, in general, different from g{a). This gd{c) will be approximated by the 
histogram of the vorticity distribution computed from a numerical simulation. 
While g{a) would be conserved by an inviscid dynamics, only gdicr) may be 
constructed from the computed field, and will vary in time. Only at the initial 
time, by definition, g{a) = gdic). A consequence of the DVC is that the 
equilibrium field ujq is the one that has maximal energy among all the fields 
with the same distribution gd{c). 



2. 6 Applications of the maximum entropy theory 



Several recent papers compare direct numerical simulations with the predic- 
tions of the maximum entropy theory. Even recognizing their value, we think 
that their Ansatze are not justifiable in our cases, or that their conclusions 
do not match our results. In detail, Miller et al. (1992) take into account a 
two-valued (0,C(Jo) and a three-valued (0,c<Jo,— c^o) initial condition onto a cir- 
cular corona with no slip boundaries. They solve the variational problem by 
a Montecarlo dynamics and compare it with a direct numerical simulation of 
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the flow. Whitakcr and Turkington (1994) consider two equal circular vortex 
patches on a closed disk, with zero ambient vorticity. They use an iterative 
solver for the constraint equations, and compare the result with more exten- 
sive contour dynamics simulations on the infinite plane (not on the disk). 
Sommeria et al. (1991) consider a shear layer, with periodical boundary con- 
ditions in the x direction and shp walls in y. Their simulation starts from a 
two-level initial condition (uj = 0,a), and the nonlinear eigenvalue problem is 
solved accordingly. In a particular limit, they are able to compute analytical 
solutions, which exhibit a bifurcation in the parameter space. One of those 
branches corresponds to solutions which break the symmetry in x, but are 
stationary (except for a system of reference mean velocity) • Their a-posteriori 
fit of the uj{ip) scatterplot is satisfactory only for one tract of the complete 
curve. They also mention simulations of multiple shear layers, saying that dif- 
ferent local vortices are obtained, preventing the system to achieve a steady 
state ("the system tends to achieve local equihbria into vortices faster than 
the global equilibrium"). Additional simulations of an idealized jet are done 
by Thess et al. (1994); also in that symmetry breaking solutions are 

found from the maximum entropy theory, and a better quantitative match is 
obtained. Symmetry breaking for periodic square, periodic channel and box 
boundary conditions is further discussed along those lines in a subsequent 
paper (Jiittner et al. 1995). 

In the maximum entropy setting, continuous symmetries generate additional 
terms in the exponentials of the equations (19). If these terms involve an ex- 
plicit dependence on the coordinates, such as in the case of the conserved 
angular momentum on the infinite plane, the soultion of (19) could be non 
stationary (Robert and Sommeria 1991). In the periodic case, however, no 
continuous symmetry besides the translation exists, and this possibility is 
prevented. 

Chavanis and Sommeria (1996) start assuming that u!{ip) is linear, as it is in the 
minimum enstrophy context, and give analytical maximum entropy solutions 
for rectangular and circular closed domains. This is said to be justifiable in 
a particular 'strong mixing' limit. The solutions are always stationary, but 
admit monopole/multipole bifurcations, which are thoroughly listed. 

A true attempt to validate the maximum entropy rather than the minimum 
enstrophy theory in an experiment is done by Huang and Driscoll (1994). 
They consider a rather simple metaequilibrium profile of a magnetized electron 
column, which obeys to the 2D Euler equation. They conclude that the mixing 
cannot be assumed to be ergodic and that the closest fit to the data is provided 
by the numerical minimum enstrophy solution. 
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3 Numerical experiments 



We performed a number of numerical experiments. To this extent we inte- 
grated in time several vorticity fields, using a rather standard protocol for the 
Navier-Stokes equation. We used a two-dimensional 2/3 dealiased pseudospec- 
tral code on the periodical square (27r)^. This choice, which fixes the Green 
function of the problem, was done for easiness of implementation. Nothing 
in the previously exposed theories prevents us to use this particular choice 
of boundary conditions. The Green function becomes G(k) = in Fourier 
space, as known. The (undealiased) resolution of the vorticity fields considered 
ranges between 16^ and 512^. A small viscosity is introduced mainly for numer- 
ical purposes, in order to prevent finite size effects, such as high wavenumber 
pile-up. Viscosity is seen to be too small if the energy, as seen from the E{k) 
spectrum, accumulates at high k. According to a generally accepted practice, 
we used hyperviscous dissipativc terms of the form U2V^uj or usV^u in place 
of the ordinary viscous term. Even if drawbacks in the use of hyperviscosity 
instead of normal viscosity have been recently discussed (Jimenez 1994), we 
think that the choice is in practice uninfluent to our simulation. We confronted 
runs in which the same initial conditions were integrated with different forms 
of the dissipative term, observing almost no influence on the resulting final 
configuration. The small scale character of the hyperviscosity suppresses also 
some truncation effects, such as Gibbs wiggles in the proximity of steep gradi- 
ents of the vorticity. Some dissipation is needed to 'underresolve smoothly' the 
smallest features which form during the evolution. The amount of dissipation 
has to be grossly matched with the rate of creation of smallest-scale structures, 
but does not have to be fine tuned. Time marching is accomplished by a fourth 
order Rungc-Kutta integrator, with explicit treatment of the dissipativc term. 
A Courant-Friedrichs-Levy condition is employed to vary the time- step. The 
adoption of At — O.SAx/IVmax] appears to be accurate enough. 

Most of our runs started from initial conditions consisting of variously ar- 
ranged constant patches of vorticity. In most cases wc simply distributed equal 
areas of vorticity equal to ±1 on the square. This already provided us with 
quite a variety of outcomes, and allows direct connection with the various for- 
mulas of section 2. A pattern of this sort is the family of the 'fuzzy checkers'. 
By these we mean checkers with randomly perturbed edges. The perturba- 
tion is done in a way that insures a constant zero mean vorticity. A regular 
checker of vorticity is an unstable stationary field (the velocity is orthogonal to 
the gradient of vorticity, which is significantly different from zero only on the 
boundaries), while 'fuzzy' checker is immediately destabilized. Another possi- 
bility to generate (usually) unstable initial conditions, is to consider random 
arrangements of rectangular ur = ±1 tiles upon the domain. A little smooth- 
ing on the initial conditions is actually necessary to prevent the edge ringing 
mentioned above, but, again, it does not seem critical for the results. As re- 
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.198 
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Table 1 

Resolution, hyperviscous coefficient fs, total integration time Tj, initial and final 
values of energy and enstrophy Eq, Ef, Qq, Qf for the runs presented. The last two 
columns give the increase in energy of the "prototype rearrangements" described in 
Appendix A. 3, relative to Ef. o stands for the 'square dipole' and = for the 'smooth 
stripe'. 

marked above, those initial conditions, even if a little odd in appearance, are 
perfectly legitimate as test cases for the relaxation models. The resolutions of 
the various runs presented, together with the hyperviscous coefficient 1/3, the 
total integration time Tf, the initial and final values of energy and enstrophy 
£■0, Ef, Qo and Qf, are reported in Table 1. 

A typical case is shown in figure 1. The patches initially intertangle in a com- 
plicated way. The process generates a lot of filamentary features, which, both 
because of the finite resolution and the small-scale dissipation, are blurred out 
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and disappear from the landscape. At later times, the final shape becomes 
apparent. The process of stretching and uniformization of the details contin- 
ues, and a 'simpler' state stabilizes. This configuration undergoes no further 
evolution, apart of a slow erosion of its contours, due to the viscosity. Since 
the viscosity can in principle be made very small, this state is long lived and 
can be named 'final'. At this stage, the viscous term of equation (1) is some 
orders of magnitude smaller than the vorticity (or the strain) field, everywhere 
on the domain. The evolution of the velocity field is therefore still dominated 
by the nonlinear terms, which balance to make the convective term null and 
the vorticity configuration stationary. The timescale for viscous damping is 
much longer than the characteristic turnover time. 

In the case of figure 1 the final state is a dipole. Something similar is also 
observed in other cases. We document a few of them in figure 2. The broad- 
est vorticity contours in the final state appear squared because of the bi- 
periodic boundary conditions. This refiects the shape of the Green function 
(Glasser 1974,Seyler, Jr. 1976). The maximum and the minimum of vorticity 
are unique, and displaced each other by half of the box size in both coordi- 
nates. The final state is stationary, and this agrees with the traditional point 
of view, but not entirely. The profiles of the positive and negative cores may 
differ. The scatterplot of (wjip) tends to a line, but does not match any of 
the expected functional relations, such as (8), (13) or (22). In all cases, the 
scatterplot of {uj,iIj) is seen to lie in the first and third quadrants. This is 
consistent with any of the proposed a;(-?/') relations, in which (3 is negative. In 
addition, the different cases mix differently the available vorticity. If we plot 
the global distribution of the vorticity (third panels of the graphs in figure 2), 
we see that the final histograms vary from case to case. The profile of the final 
vortex cores is therefore peculiar to each decay. In all of them, however the 
final vorticity has a single maximum and sigle minimum, and the dipole may 
be thought as an arrangement of the available vorticity around two centers. 

The final dipole is not the only possibility, though. In figure 3, we see a differ- 
ent outcome. In this pattern of two stripes is formed. The boundaries 
of these two stripes translate in opposite directions varying slightly their cur- 
vature according to their relative position. The state is therefore recurring. 
Interestingly, the two zones cannot be explained simply by the absence of 
mixing of the original patches: detailed analysis shows that each stripe con- 
tains a fraction of the fiuid of originally opposite vorticity. The mixing inside 
a single stripe is complete, so that the vorticity is almost constant and lower 
than that of the likely-colored patch in the initial condition. The effectiveness 
of mixing within the stripe, but its absence among the two stripes, is clearly 
evidenced when including numerical passive markers into the flow (flgures 
omitted) . 

In figure 4 we see other examples of fields which decay in nonstationary final 
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states. There is a variety of forms: non-symmetric fat cores, which oscillate 
in time; wavy stripes with central blobs that translate in one direction. Yet 
the last case shown in figure 2 is not really a simple dipole; actually three 
smaller and long lived vortices happen to form inside the main core. They 
keep orbiting for long time and weaken very slowly because of the viscosity, 
but arc never able to merge. It is sometimes common, as in the case of figure 
3, to get a complete uniformization of the vorticity in different well-separated 
zones of the flow. A number of patches of almost constant vorticity are thus 
formed, and their arrangement prevents further mixing. In all these cases the 
scatterplot of never thins out, and, a fortiori cannot be fitted by any 

of the relations given in section 2. Nevertheless, these configurations are long 
lived. Several features are still common to both types of decays. The final 
states are "dipoles" in a broad sense, in that they show two simply connected 
regions of vorticity above and below the mean. The local maxima of |ci;| are 
also absolute maxima. The small-scale structures are generated (and subse- 
quently blurred) quite early in the simulation. If we compute characteristic 
eddy turnover times based on the sizes of the patches in the initial conditions, 
we see that the palinstrophy (associated to the amount of details) reaches 
its maximum within a few turnarounds. The entire simulation is carried on 
typically for some tens of turnover times. The loss of energy during the decay 
is always negligible. The enstrophy is seen to decay significantly during the 
decay, until the final state is formed; it then remains nearly constant. Most of 
the theories proposed in section 2 take into account a lowering of Q, but fix it 
in an unique way. We see that not even the initial E and Q are sufficient to 
determine the final state: as a counterexample, the case of figure 3 and the one 
shown in the second row of figure 4 possess initially nearly the same values 
of E and Q, and the same distribution of vorticity, but reach different final 
states. Other counterexamples may be easily generated taking an intermedi- 
ate configuration and reversing the time integration. The blurring of details 
is an irreversible process, and the backward integration does not lead back 
to the initial condition. The forward and backward final states have generally 
different properties. 

For the initial conditions we have used, the final states appear to be quite inde- 
pendent from the resolution and the (small) values and forms of the viscosity. 
Generally, a higher spatial resolution just requires a longer time to finish to 
smooth the details and to reach the final state. To ground this affirmation we 
give some examples. 

First, we repeated the simulation shown in figure 3 with four different values 
of the hyperviscosity. The appearance of the field after long times is shown 

in figure 5. The same wavy pattern, with uniform stripes of vorticity and 
undulated boundaries is seen. Phase differences between the final panels are 
not meaningful, since the field is plotted at different times. To demonstrate 
the convergence and to provide one example of the evolution in time of the 
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characteristic quantities, we plot, in figure 6, E[t), Q{t) and P{t) for tlie 
different values of z/3. The functions E and Q always stabilize on the final 
level. The loss in E is indeed negligible, less than 1% in the most dissipated 
case. The enstrophy instead decreases of roughly a factor two, but stabihzes 
to a final value independent of the hyperviscosity. The palinstrophy is seen to 
increase in time in the initial stages. This behavior, which is also predicted 
by closure theories' estimates, is related to the increase of small-scale details. 
The maximum value is achieved when such details begin to be more rapidly 
numerically underresolved than generated. This can be seen by looking at the 
field in physical space (compare for instance the intermediate panels of figure 
3). The curves are also seen to scale inversely to the hyperviscosity. This 
behavior, apart from different values and different relaxation times, is also 
found in all the other decays, both the stationary and nonstationary ones. It 
is also observed when normal viscosity or second viscosity are employed. 

Secondly, in order to show that the final state found is not an artifact of the 
finite resolution, we display the integration of another arbitrary initial condi- 
tion at different spatial resolutions (figure 7). The hyperviscosity coefficients 
are adjusted case by case according to the highest wavenumber retained. We 
see that unless the resolution is extremely poor (16^), a final state with given 
characteristics is obtained. In this case, the final state is nonstationary, the 
two elongated cores translate in one direction, while both background stripes 
convect in the opposite one. 

We note, passing by, that the Fourier spectrum of these final states is always 
decreasing in k. We are not however concerned with possible laws to fit their 
slopes, and we do not report them here. In general, the final state is neither 
a pure |k| = 1 state (as predicted by the minimum enstrophy principle), nor 
a spectrum introduced in section 2.1. Most of the energy of the final field is 
contained in the first modes, but this fact itself does not explain the formation 
of a huge variety of profiles. 



4 Conclusions and perspectives 

We have considered several examples of relaxations of two-dimensional vortic- 
ity fields, and compared them with the theories which have been so far pro- 
posed in literature. We have brought some numerical evidence that contrasts 
with the existing final states theories. The different theories do not agree in 
their predictions, as can be seen in the different a;o('0o) which they propose. We 
have shown the existence of legitimate final states which are not included in the 
accepted point of view. Nonstationary ones are among them. The remark that 
stable, nonstationary configurations may indeed form from unstable ones is 
indeed not new. Stable tri- and quadrupole satellite systems, which result from 
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the decay of unstable vortex blobs, have already been investigated analytically 

and numerically (Carton et al. 1989, Carton 1992, Morel and Carton 1994, Carton and Legras 1994). 
Beautiful experimental evidence is also provided (van Heijst et al. 1991, Flor and van Heijst 1996). 
In all these cases {uJ,ip) is derived from the experiment, and is seen to evolve 
from a non-monotonous to a branched monotonous curve. On the other hand, 
also the non- uniqueness of stable dipolar solutions on the periodic box has been 
stressed (Rasmussen et al. 1996). Numerical experiments shown there display 
{uj,ip) scatterplots which are tentatively fitted by polynomial relations. For 
this reason we have not tried to fit one or the other of the free parameters of 
the reported models versus our numerical results: the possibility of a good fit 
appears more accidental than general. In our opinion many of the underlying 
assumptions of these theories are not always granted, and call therefore for 
a more careful treatment. For instance, we think that the final state cannot 
be predicted by a 'final state theory' which ignores completely the dynamical 
path underwent by the relaxing system. This is demonstrated by the fact that 
configurations with the same two initial vorticity levels decay to different final 
states. The quest for appropriate statistical descriptors, which are preserved 
during the evolution, is still open. For instance, it is not known how the initial 
macrovorticity distribution relaxes to the final one. It would be appealing to 
parametrize this change by a few dynamical quantities, such as the decrease 
of Q. We have supposed that some of the shapes of the final states can be 
understood as rearrangements of the vorticity which is available at late times. 
This seems to apply at least to the final stationary states, for which the re- 
quirement of energy maximization appears as sufficient. The instability of the 
intermediate vorticity configurations is always such that the system is led to 
mix its vorticity and drawn toward the final state. In the other cases we may 
conjecture that, for some stability reason, the available vorticity is unable to 
collapse onto a stationary configuration. In many of the previous works, initial 
Gaussian fields with random spectra were considered. These fields show gen- 
erally a smaller population of absolute higher vorticity (i.e., Qdic) with long 
tails). It is often observed (but not always) that a lot of concentrated vorticity 
patches are formed in the early times of evolution and they survive relatively 
long times (see e.g. Brachet et al. 1988, Kida 1985). Such intermediate states, 
composed of several isolated cores, which have so often been described in lit- 
erature, could be understood as metastable preliminary local rearrangements 
of the vorticity. 

In other cases, however, the attainment of stationary states seems prevented. 
Possibly, the system is trapped in a metastable state, which is not stationary 
(in the cases encountered in this work, periodical in time). The stationary 
state could be missed just because the dynamics of the vorticity is by itself 
insufficient to trigger a catastrophic mixing process, which alone could alter 
further the vorticity distribution, and lead the system to its very final state. It 
would then be very interesting to find the requirements for this met ast ability. 
An implication of such metastabilities would be that some kind of external 
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forcing might trigger the transition to the stabler states and some perhaps 
not. 
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Appendix: Energy nicLximization and stability 

All the theories proposed in section 2 identify the final state as the one which 
solves a certain extremum problem. In the case of the minimum enstrophy 
principle, one may think to rescale the vorticity, and to rephrase the problem as 
to the search for the maximal energy, given a fixed enstrophy. In the setting of 
the minimum entropy, the DVC instead tells us that the final state is actually 
also a state of maximal energy (with fixed ga). It is therefore worth to add 
some remarks on energy maximization. 

We are not able to determine in how many equally energetic ways the same 
amount vorticity can be distributed over a domain. We might ask if equally 
energetic states which are not dynamically related, i.e., which do not evolve 
one into the other, because they are topologically different, could converge 
to the same final state. We cannot answer to those questions, but we can 
add a few considerations which at least rule out some possibihties. These 
provide necessary requirements, but they do not determine the configuration 
of vorticity nor the fiow. 

A.l Infinitesimal deformations 

Let us first consider a generic area preserving infinitesimal deformation on a 
two-dimensional domain. The deformation may be written as Sx — V"'"50(x), 
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where is an arbitrary field. It is easily seen that such transformations preserve 
all the moments of the vorticity: 



5 J JdF^^l j J{u,5(t))J-^(f 

^-ijdcp [K^'-^)^ - {ujyUj'-%] ci'x = 0, (28) 

Integration by parts is carried out since the boundary terms vanish for any 
boundary condition. The variation of the energy with respect to such defor- 
mations is 

SE^S^ ya;(x)a;(x')G'(x',x)dWx' = y 5u;(xV(x')G'(x', x)dWx' = 
= - y J (50(x), cu(x)) c^(x')G(x', x)dWx' = 

= y 50(x) J (u;(x), ^(x', x)) a;(x')d^xdV . (29) 
Therefore the extremal energy states satisfy 

y J (a;(x), ^(x', x)) a;(x')(iV = J (a;(x), ^(x)) = 0. (30) 



5E 
50 



Thus if the dissipation is absent, the stationary states extremize the energy 
with respect to infinitesimal incompressible deformations. These states are 
"local" extrema. 



A. 2 Pointwise mixing exchanges 



We can write a transformation which transfers some of the vorticity present in 
the neighborhood of Xi to a similar neighborhood of X2 and vice versa. To this 
extent we employ a shape function 5e(x), which is equal to 1 in a neighborhood 
of radius e of x = (modulo the boundary conditions) and zero everywhere 
else. The vorticity field is transformed according to 



a;(x) a;(x) + r [c^(x - xi + X2) - c^(x)] (5e(x - Xi) 

+ r [a;(x + xi - X2) - u;(x)] 6^{x - X2) . (31) 

The parameter r may vary from to 1; in the latter case, an exchange between 
the vorticity in Xi with that in X2 is realized. The energy of the transformed 
field is computed by substituting the above expression in equation (3). When 
consider the limit e ^ 0, the integrals are approximated by the mean value 
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of the integrand times the area of the integration region, provided that u 
is continuous. The terms multiphed by r come of order O(e^), while the ones 
with remain 0{e^). The singularity in the Green function poses no problems, 
since for small arguments G(x, y) ~ — In |x— y|, which is integrable. Therefore 

E ^ £; + r7re'[V'(x2) - V(xi)][a;(xi) -a;(x2)] +r'0(e^) . (32) 

This gives us the infinitesimal change in energy due to a transformation which 
modifies infinitesimally, but not continuously, the field. The fact that this 
transformation is not dynamical is not of our concern, as far as we look for 
properties of the configurations and not for their evolution. 

If we refer to the maximal energy state, its energy shall decrease whichever 
the couple of points Xi and X2 and the value of r. Therefore, cj(x2) > cj(xi) 
implies ^(x2) > V'(xi) and vice versa. The most important consequence of it is 
that ip{uj) must be single valued and monotonous. This is seen by considering 
the plot of a;(x) versus The monotonicity is imphed by the fact that 

data points on the (a;, ■0) plane can be labelled by x. If the relation uj{;ip) was 
non monotonous or more than single valued, then it would be possible to find 
points at which a;(x) < a;(xi) and '^(x) > '^(xi), contrary to the hypothesis. 
As a further consequence, the absolute maximal energy state is also stationary, 
because a;(x) and are functionally related. 

A particular case is achieved when the points are infinitesimally close and 
displaced in the direction d: in this case (32) is read as (d • V'?/')(d • Va;) > 
at any point and along any direction. 

Related to this, but derived with a different machinery, is the Rayleigh- Arnold 
criterion (Holm et al. 1985,McIntyre and Shepherd 1987,Carnevale and Vallis 1990). 
It states that among all the isovortical fields the maximal energy arrange- 
ment is unique and stationary. Furthermore, if dipo{u!o) / dcuo, or equivalently 
—Vipo ■ Va;o/| Va;oP are strictly positive and limited, then the field is also non- 
linearly stable. We also note that there may exist also other stable states, in 
particular those obtained by means of some symmetry transformation, which 
is specified by invariants inexpressible as functionals of the vorticity. Ana- 
lytical stabihty criteria other than Arnold's one, that refer to local features 
in ambient fiows (e.g. Nycander 1995), do not appear to be of direct use for 
assessing global stability. 

A. 3 Energy maximization algorithms 

We have seen in A.l and A. 2 properties of 'energy-maximizing arrangements 
of the available vorticity'. By 'arrangement', we mean another field, which 
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has exactly the same fractional area for each level of vorticity (i.e. the same 
vorticity histogram), but a different configuration. The preceding facts tell 
us that these rearrangements are unique, stationary and stable. We are not, 
however, guaranteed that any initial field relaxes into such an energy maxi- 
mizing configuration. Not only it is a priori undetermined how the histogram 
of macrovorticity evolves in time, but also not all fields are seen to approach 
a stationary final state. We want therefore to characterize the maximal en- 
ergy state for a given vorticity population, and its relation with the observed 
dynamical final state. 

One possibility would be provided by the Carnevale, Vallis and Young pseu- 
dodynamics (Vallis et al. 1989, Carnevale and Vallis 1990). They proposed a 
numerical procedure to deform a vorticity field by increasing its energy and 
conserving its enstrophy. The procedure does not explain the geometrical prop- 
erties of the extremal states, which can only be said to be local, and not ab- 
solute maxima of the energy. Lacking analytical procedures, we rather used a 
simpler alternative combinatorial approach. Starting from a given configura- 
tion of vorticity, two square tiles are randomly exchanged and the new energy 
is computed. The new configuration is kept if the energy increases, discarded 
otherwise, and the process repeated. Miller et al. (1992) also use a more sophis- 
ticate numerical procedure of this kind. Such approach is easy to implement, 
but very slow in convergence when the number of tiles is significant. With the 
aid of this tool, we observed that the highest energy configuration achieved 
by permutation of equal area tiles of vorticity ±Ui on the periodic square, 
is a subdivision in two stripes parallel to the sides. When intermediate levels 
of vorticity are present, as in our runs, where patches' edges are significantly 
smoothed, the most energetical of these permutations of tiles has sometime 
still a striped shape, sometimes a central cored appearance. 

To cope with high resolution rearrangements, we chose to limit ourselves to 
compare energies of only a few selected prototypical configurations. To pro- 
duce them we sorted the discrete tiles of the field according to value of their 
vorticity, and deployed them upon the domain following the same order of a 
sample configuration. In this way we obtain a rearrangement with the same 
set of contour lines of a given prototype. We can then check which, among few 
prototype rearrangements, increases the current energy. One of the test profiles 
we chose is the "square dipole", which is modeled on G'(x) — G (x + (vr, vr)); its 
contour lines look similar to those of various dipolar final states shown in Fig. 
2. Another one is the parallel stripe, which is the maximal arrangement for a 
two-level population. If we compute the energies of the rearrangements of the 
final states obtained in section 3, we basically see that for the stationary final 
states the highest energy is achieved for the 'square dipole' rearrangement, 
with minimal increase of energy. For the nonstationary final states, such as 
the wavy ones, often the stripe rearrangement results more energetical. The 
actual values of the energies are reported in Table 1. The choice of the max- 
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imal profile is basically a function of the available vorticity. Still, the energy 
of the nonstationary states is often significantly increased by the procedure. 
This indicates that the final state, in those cases, is far from the highest energy 
configuration for the given g^. Conversely, the final state may be dynamically 
inaccessible from the rearrangement: such is the case for the wavy configu- 
rations, if we maintain that the energy is almost conserved. In addition, the 
wavy cannot be seen as a perturbed stripe, since the maximal energy state 
cannot be unstable to transversal perturbations, thanks to the Arnold stability 
criterion. 

We are not able to provide a full analytical treatment of the family of the recur- 
ring wavy solutions here, but we would like to add a small piece of numerical 
evidence. We plot the energy corresponding to all two level striped fields with 
sinusoidal undulations, with varying amplitudes and phase shift (figure 8). We 
see that the energy has to decrease significantly with increasing amplitude in 
order to allow for undulations of the boundaries. For a range of energies below 
the highest, energy isolines span the entire interval < A(f) < 27i, with the 
amplitude depending on the phase shift. A recurring evolution of the field, 
with pulsating amplitudes, is therefore allowed. Without attempting to make 
a precise fit to our final states, we claim that the picture is appropriate. It is 
clear from the figure that for amplitudes of the undulation higher than A ~ jTT, 
but still lower than the geometrical limit A = ^tt, the range < A0 < 2tt is 
no more entirely accessible at constant E (the separatrix is drawn in bold) . In 
fact, numerical tests with initial conditions of this kind exhibit stability below 
an amplitude threshold, and a catastrophic mixing above. 
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Fig. 1. A typical case: decay of a two-level patched vorticity field into a stationary 
dipole. Plots of the vorticity field at various times. A few vorticity isolines are drawn 
for reference. 



Fig. 2. Arbitrary 2-level initial conditions which decay in final dipoles. Each row 
illustrates one case, with, in each panel: the initial and the late time vorticity field; 
the initial (dashed line) and final (full line) global vorticity distributions; the final 
scatterplot of {tp,oj)- All cases start from initial patches with a; = ±1, except for the 
last one, where u = —1.5, +0.5. The mean vorticity is always zero. 



Fig. 3. A less typical case: decay of a two-level patched vorticity field into a wavy 
pattern. 



Fig. 4. Arbitrary 2-level initial conditions which decay in nonstationary configura- 
tions. The format of presentation is the same as in Fig. 2. 



Fig. 5. Decay to the final wavy state, obtained with different hyperviscosities. The 
format of presentation is the same as in Fig. 2. The value of is respectively 2 -10"^, 
2 • 10"^, 2 • 10"^°, and 2 • 10"^^ 



Fig. 6. Energy, enstrophy and palinstrophy decay for different values of the hyper- 
viscosity. 



Fig. 7. Effect of lowering the resolution of the simulation. The same initial condition 
is integrated at 512^, 256^, 64^, 32^, 16^, and hyperviscosity rescaled appropriately. 
The format of presentation is the same as in Fig. 2. 



Fig. 8. Energy of a family of wavy patterns of vorticity equal to ±1, as a function 
of A and A0. The region of the (A, A(j)) plane where the two stripes would interfere 
is drawn cross-hatched. 
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